%% 2D Ising Model

%web('http://en.wikipedia.org/wiki/Square-lattice_Ising_model')

clc
close all
clear all

Beta=1;         %1/KbT
J=-.2;            %Interaction Strength
B=0;            %Magnetic Field 
N=100;           %Number of spins

%There is phase transition at -1/(Beta*J)=2/log(1+sqrt(2))

-1/(Beta*J)

%Start with random orientation
Spins=double(randn(N,N)>0);
Spins(Spins==0)=-1;

%Plot Spins
FigH=figure;
%FigH.Position=[200,200,900,400];
imshow(Spins,[],'InitialMagnification',400)
s=sprintf('Beta=%g, J=%g, B=%g',Beta,J,B);
title(s)

%Setup so we stop when figure is closed
FigH.CloseRequestFcn='Go=0;delete(FigH)';
Go=1;

%Run Markov Chain Monte Carlo
NChain=1000000;
nn=0;
while Go
    nn=nn+1;    %counter for display 
    
    %Make test jump
    
    %Pick site and flip spin
    IDX=ceil(N*rand);
    IDY=ceil(N*rand);
    SpinsTest=Spins;
    SpinsTest(IDX,IDY)=-Spins(IDX,IDY);
    
    %Energy from nearest neighbor interactions
    E_Spin=J*(sum(sum(Spins(1:end-1,:).*Spins(2:end,:)))+... %Y neighbor
        sum(Spins(1,:).*Spins(end,:))+... %Y periodic BC 
        sum(sum(Spins(:,1:end-1).*Spins(:,2:end)))+... %X neighbor
        sum(Spins(:,1).*Spins(:,end))); %X periodic BC
    
    E_SpinTest=J*(sum(sum(SpinsTest(1:end-1,:).*SpinsTest(2:end,:)))+... %Y neighbor
        sum(SpinsTest(1,:).*SpinsTest(end,:))+... %Y periodic BC 
        sum(sum(SpinsTest(:,1:end-1).*SpinsTest(:,2:end)))+... %X neighbor
        sum(SpinsTest(:,1).*SpinsTest(:,end))); %X periodic BC

    %Energy from Interactions with feild
    E_Field=-B*sum(Spins(:));
    E_FieldTest=-B*sum(SpinsTest(:));
    
    %Total Energy (Hamiltonian)
    H=E_Spin+E_Field;
    HTest=E_SpinTest+E_FieldTest;
    
    %Test acceptance probability for Metropolis-Hastings
    Alpha=exp(-Beta*(HTest-H));
    
    %Metropolis-Hastings rule
    if rand()< (Alpha)
        Spins(IDX,IDY)=-Spins(IDX,IDY);
    end
    
    ShowIter=1000;
    %Show every ShowIter interations
    if nn/ShowIter==round(nn/ShowIter)
        figure(FigH);
        imshow(Spins,[],'InitialMagnification',400)
        pause(.01)
    end
end





